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ABSTRACT 

To assess the effect of baryonic "pinching" of galaxy cluster dark matter (DM) haloes, 
cosmological (ACDM) TreeSPH simulations of the formation and evolution of two 
galaxy clusters have been performed, with and without baryons included. 

The simulations with baryons invoke star formation, chemical evolution with non- 
instantaneous recycling, metallicity dependent radiative cooling, strong star-burst, 
driven galactic super-winds and the effects of a meta-galactic UV field, including 
simplified radiative transfer. The two clusters have T~ 3 and 6 keV, respectively, and, 
at z both host a prominent, central cD galaxy. 

Comparing the simulations without and with baryons, it is found for the latter 
that the inner DM density profiles, r < 50-100 kpc, steepen considerably: Aa ~ 0.5- 
0.6, where -a is the logarithmic DM density gradient. This is primarily due to the 
central stellar cDs becoming very massive, as a consequence of the onset of late time 
cooling flows and related star formation. Once these spurious cooling flows have been 
corrected for, and the cluster gravitational potentials dynamically adjusted, much 
smaller pinching effects are found: Aa ~ 0.1. Including the effects of baryonic pinching, 
central slopes of a ~ 1.0 and 1.2 are found for the DM in the two clusters, interestingly 
close to recent observational findings for galaxy cluster Abell 1703 based on strong 
gravitational lensing. 

For the simulations with baryons, the inner density profile of DM and cluster gas 
(ICM) combined is found to be only very marginally steeper than that of the DM, 
Aa < 0.05. However, the total matter inner density profiles are found to be Aa ~ 0.5 
steeper than the inner profiles in the dark matter only simulations. 

Key words: cosmology: theory — cosmology: numerical simulations — galaxies: 
clusters — galaxies: formation — galaxies: evolution 



1 INTRODUCTION 

Large N-body cosmological simulations have been carried 
out for a decade, with the goal of making statistical pre- 
dictions on dark matter (DM) halo properties. Because of 
numerical issues, most of these large cosmological simula- 
tions contain dark matter particles only. They all reliably 
predict that the 3D density profile pdm('') should fall as 
at large radii (but within the virial radius). Obser- 
vations have confirmed these predictions for cluster sized 



vations nave connrmed tnese predictions tor cluster sized 
haloes (e.g. [K ncib ct al.' '2003"; 'Pointcco uteau et al.l |2005| ; 
iMandelbaum et al .. .2008, ; .Okabc ct al.200d ) . This agreement 
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is likely to be connected with the fact that at large radius, 
the density profile of a galaxy cluster is dark matter dom- 
inated and the influence of baryons can be neglected. On 
smaller scales, parameterizing the 3D density profile of the 
DM using a cuspy profile pdm oc r~"; dark matter only 
simulations predict a logarithmic slope a ~ 1-1.5 for r — > 0. 
The ex act value of the cent r al slope and its universality is de- 
bated jNavarro et al.lll997l;lMoore et al.lll998l:lGhigna et all 



I2OOOI ; iRicottil I2OO3I ; iNavarro et all 12004 ; ICao et all 120081 ) 

Theoretical efforts indicate that the inn er logarithmic slope 
of DM halos should be a ~ 0.8 (|Austin et all 120051 ; 
[Hansen fc Stadeill2006l ). Although clearly of interest, these 
dark matter only studies and their predictions do not help 
much to make comparison with observations. Indeed, ob- 
serving the central part (i.e. the inner ~ 500 kpc) of a galaxy 
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cluster at any wavelength reveals the presence of baryons 
(in the forms of stars and X-ray emitting hot gas). Thus 
any attempt to compare observations to simulations in the 
center of galaxy clusters has to be made with numerical sim- 
ulations (or calculations) taking into account the baryonic 
component and its associated physics. 

On the numerical side, efforts are currently being made 
to include effects of baryons in a realistic way in the simu- 
lations (see below). Problems remain, however: for example, 
due to the so-called "over-cooling problem" , numerical sim- 
ulations tend to predict central brightest cluster galaxies 
which are too blue and too massive relative to observations. 
Different effects do compete when it comes to the central 
slope of the density profile: the cooling of gas in the central 
regions of galaxies and clusters is expected to lead to a more 
concentrated dark matter distr ibution (th e so-cal led adia- 
batic contraction, see ^ l uni cnt hal et aDll98 6: Gn cdin et al.l 



I2OO4I : iGustafsson et alTl2006[ '). On the other hand, dynam- 



ical friction heating of massive galaxies against the diffuse 
cluster dark matter could in principle fla, t ten the slope of the 
DM density p rofile (El-Zant ct al.l l200fl : iNipoti et al.1 [20031 . 
I2GO4I : iMa fc B ovlan-Kolcliin .200'l). Also, the properties of 



the inner part of simulated galaxy clusters (even in dark 
matter only simulations) can de pend significantly on in i- 
tial conditions as demonstrated in lKazantzidis et al.l (|2004 ) . 
Consequently, no coherent picture has yet emerged from N- 
body simulations when it comes to the shape of the inner 
density profile of structures. 

On the observational side, efforts have been put on 
probing the central slope a of the underlying dark matter 
distribution. These analyzes have led t o wide-ranging re- 
sults, whatever the m e thod used: X-ray (lEttori et al. 20021: 
Arabadiis et al.1 12OO2I: iLewis et al.l 120031; IZappacosta et a ] 



20061): lens i ng (iTvson et al.1 1 19981 : ISmith et al.1 
Dahle et al.1 l2003l: ISand et al.1 I2OO2I: ICavazzi et al 



Gavazzil |200"5[ 



Limousin et al 



Sand et al.l |200. 



Bradac ct al 




20081: iRichard et alii2009l: lOeuri et at 120091 ) 



or dvnamics (Kelson et al.ll2002l : iBiviano fc Saluccill200^ . 
This highlight the difficulty of such studies and the possible 
large scatter in the value of a from one cluster to another. 

In summary, one needs to probe observationally and nu- 
merically the behavior of the underlying dark matter distri- 
bution (i.e. after the baryonic component has been separated 
from the dark matter component) in the central parts of 
galaxy clusters. The main difficulties are: i) Observationally, 
to be able to disentangle the baryonic component and the 
underlying dark matter distribution; ii) Numerically, to im- 
plement the baryonic physics into the simulations; iii) Then 
to c ompare bot h approaches in a consistent way. 

ISand et al.1 (|200l |2004| . l2008l ) carried out lensing ana- 
lyzes aiming to probe the central mass distribution in six 
galaxy clusters, using the measured velocity dispersion pro- 
file of the cD galaxy as an extra constrain. They found cen- 
tral density slopes smaller than 1. 

Recently, there have been a growing interest on 
AbeU 1703, a massive 2 = 0.28 (|Allen et al.l Il992l ) X- 
ray cluster with a luminosity Lx = 8.7x10** ergs~^ 
l|Bohringer et al.ll2000l ). It contains a large number of grav- 
itational arcs, enabling a very detailed lensing analysis. 
Moreover, although it displays an intriguing filamentary 
structure along the north-south direction, it looks rather 
circular from the geometrical configuration of its multiply 



imaged systems, in particular its giant arc, located at large 
angular separation (~ 35" ~ 147 kpc). It is a uni-modal 
cluster, likely to be relaxed and characterized by a central 
dominant elliptical cD galaxy. This makes it much easier 
to interpret the results of the mo deling compared to bi- 
modal clusters such as Abe ll 1689 (iLimousin et al.ll2007bl: 
Rieiner-S0rensen et al.ll2009l), A b eU 22 18 (|Eh'asd6ttir et al.1 
20071 ). AbeU 68 (|Richard et al.1 |2007| ) or MS 2053.7-0449 
(Verdugo et al., 2007). In fact, regu l ar rel axed clusters are 
rare at such redshifts (|Smith et al.l |2005| ) . Finally, it dis- 
plays a remarkable lensing configuration, forming a cen- 
tral "ring" composed of four bright images that represent 
a ra re example of lensing by a hyperb olic umbilic catastro- 
phe (jOrban de Xivrv fc Marshal]||2009l ) . This ring is located 
close to the central cD galaxy (4-9 ", corresponding to 17- 
38 kpc), providing a robust constraint in the very central 
part of the cluster. 

This motivated several lensing works on Abell 1703 
aimed to probe the dark matter density slope of the central 
region (~ 20 — 200 kpc), modeling the dark matter distribu- 
tion of Abell 1703 by a generalized NFW profile, viz. 



p{r) 
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(1) 



where is the scale radius. These works (ILimousin et all 
l2008l : [Saha et al ] |2009l : iRichard et aill2009l : lOguri et alll2009D 

determined the inner slope aNFW to be consistent with the 
universal profile (aNFW = 1)- The m ore recent strong lens- 
ing analysis by IRichard et all (|2009l ) , based on the identi- 
fication of 16 multiply imaged systems, of which 10 are 
spectroscopically confirmed, found aivFW=0.92±0.04 (la 
confidence le v el). N ote that the scale radius derived by 
IRichard et"ail (|2009|) is found in agreement with that found 
by Oguri et all ( 20091 ) when combining strong and weak lens- 
ing. Since the scale radius needs weak lensing data to be 
properly constrained, this is important given the degenera- 
cies arising between the scale radius and the inner slope. 

Assuming that pure dark matter haloes are well de- 
scribed by the standard NFW profile with aNFW=^, the 
indication of the above results on Abell 1703 is that the 
baryonic pinching of the dark matter halo, probably mainly 
caused by the central dominant cD galaxy, does not sub- 
stantially affect the dark matter density slope in the central 
region of the cluster. 

On the theoretical/numerical front it has only re- 
cently been possible to carry out fully cosmological gas- 
dynamical/N-body simulations of the formation and evo- 
lution of galaxy clusters at a sufficient level of numeri- 
cal resolution and physical sophistication that the cool-out 
of gas, star-formation, chemical evolution and gas infiows 
and outfiows related to individual cluster galax ies can be 
modeled to a re asonable degree of realism fe.g.. Kav et al 



2004 



2OO2I: iThomas et al.. 2002: Valdarnini 2003; Tornatore et al 
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2005; D'Onghia et al 



2005 



20051. 12OO6I: ISaro et al...20 06; Mua nwong et al l 



Romeo et al.l _ _ _^ 

20061: iTornatore et al.ll2007l; iRomeoet al. 2008). 



iRomeo et al] |2005l), [Rmtico ct al. (^2001 



iD'Onghia et al.l (120051 ) and ISommcr-Larscn ct al. (200. 



presented fully cosmological simulations of galaxy groups 
and clusters. The TreeSPH code used was building on 
the code used for simulating g alaxy formation (e.g., 
ISommer-Larsen. Gotz fc Portinaril |2003| ). improved to 
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include modeling of non- i nstanta neous chemical evolution 
(|Lia. Portinari fc Carrarol l2002ah . metallicity-dependent, 
atomic radiative cooling, strong supernova, and (optionally) 
AGN, driven galactic winds and thermal conduction. 
The two clusters simulated have z—0 virial masses 
Myir ~ 3 X 10^"' and 1.2 x 10^^ Mq, one approximately 
the size of the Virgo cluster and the other of the Coma 
cluster. They were both selected to be fairly relaxed, and 
both display central prominent cD galaxies at z—0 as well 
as at 2=0.28. In this paper we reanalyze the simulations of 
the two clusters with emphasis on the dark matter density 
profiles. We also present re-simulations of the clusters 
aimed at correcting for the effects of spurious, late time 
cooling-flows. Finally, we present results of new simulations 
of the same two clusters, identical to the previous ones 
except that they were carried out with dark matter only. We 
then compare the results of simulations with and without 
baryons, to quantify the effects of baryons (especially the 
cDs) on galaxy cluster dark matter profiles. 

The paper is organized as follows: the code and the 
simulations are described in section 2, the results obtained 
are presented in section 3 and discussed in section 4, and, 
finally, section 5 constitutes the conclusion. 

Throughout the paper a ACDM cosmology with S7m ~ 
0.3, = 0.7 and a Hubble constant Ho ~ 70 kms^^ 
Mpc^^ is assumed. 



2 THE CODE AND SIMULATIONS 

The code used for the simulations is a significantly 
improved version of the TreeSPH code we have 
used previously for galaxy fo rmati on simulations 
(jSommer-Larsen. Gotz fc Portinaril 2003ll : Full details 



on the code are given in iRomeo et al.1 (|2006l ): here we recall 
the main improvements over the previous version. (1) In 
lower resolution regions (which will always be present in 
cosmological CDM simulations) an improvement in the 
numerical accuracy of the integration of the basic equations 
is obtained by solving the entropy equation rather than 
the thermal energy equation — we have adopted the 
"c onservative" entropy equatio n solving scheme suggested 
bv ISprineel fc Hernouistl (j2002h . (2) Cold high-density gas 
is turned into stars in a probabilist i c way as described 
in ISommer-Larsen. Gotz fc Portinaril (|2003l l. In a star- 
formation event an SPH particle is converted fully into 
a star particle. Non-instantaneous recycling of gas and 
heavy elements is described through probabilistic "decay" 
of star particles back t o SPH particles as discussed by 
iLia. Portinari fc Carrarol (|2002ah . In a decay event a star 
particle is converted fully into an SPH particle. (3) Non- 
instantaneous chemical evolution tracing 10 elements (H, 
He, C, N, O, Mg, Si, S, Ca and Fe) has been incorporated 
in the code following Lia et al. (2002a,b); the algorithm in- 
cludes supernova of type II and type la, and mass loss from 
stars of all masses. Metal diffusion in Lia et al. was included 
with a diffusion coefficient k derived from models of the 
expansion of individual supernova remnants. A much more 
important effect in the present simulations, however, is the 
redistribution of metals (and gas) by means of star-burst 
driven "galactic super-winds" (see point 5). This is handled 
self-consistently by the code, so we set k—0 in the present 



simulations. (4) Atomic radiative cooling depending both 
on the metal abundance of the gas and on the meta-galactic 
UV field, modeled after Haardt fc Madau (1996) is invoked. 
We also include a simplified treatment of radiative transfer, 
by switching off the UV field where the gas becomes 
optically thick to Lyman limit photons on scales of ~ 1 kpc. 
(5) Star-burst driven, galactic super-winds are incorporated 
in the simulations. This is required to expel metals from 
the galaxies and get the abundance of the ICM to the 
observed level of about 1/3 solar in iron. A burst of star 
fo rmation is modeled in the same way as th e "early bursts" 
of ISommer-Larsen. Gotz fc Portinaril (|2003l l. i.e. by halting 
cooling in the surrounding gas particles, to mimic the initial 
heating and subsequent adiabatic expansion phase of the 
super-shell powered by the star-burst; this scheme ensures 
effective energy coupling and feedback between the bursting 
star particle and the surrounding gas. The strength of the 
super-winds is modeled through a free parameter /wind 
which determines how large a fraction of the new-born stars 
partake in such bursting, super-wind driving star formation. 
We find that in order to get an iron abundance in the ICM 
comparable to observations, /wind ^ 0.5 and at the same 
time a fairly top-heavy Initial Mass Function (IMF) has to 
be used. (6) Thermal conduction was im plemented in the 
code following IClearv fc MonaghanI (|l999l ). 

In this paper we present results for two simulated 
clusters, "Virgo" and "Coma". Both systems were chosen 
to be fairly relaxed (no >1:2 merging at z <1). Virial 
masses at z=0 are approximately 2.0x10^* and 8.7x10^* 
H^^Mq and X-ray emission weighted temperatures 3.0 and 
6.0 keV, respectively. The clusters were selected from a cos- 
mological, DM-only simulation of a fiat ACDM model, with 
fiiv/=0.3, 176=0.036, h=0.7 and (78=0.9 and a box-length 
of 150 h~^Mpc. Mass and force resolution was increased 
in, and gas particles added to, Lagrangian regions enclos- 
ing the clusters. The Virgo cluster was then re-simulated 
using 2.3 million baryonic-|-DM particles with mgas=m,= 
3.1x10^ and mDM=2.3xlO* h~^MQ for the high resolution 
gas, star and dark matter particles. Gravitational (spline) 
softening lengths of egas=e* = 1.4 and eDM=2.7 /i~^kpc, re- 
spectively, were adopted. The more massive Coma clus- 
ter was re-simulated using 1.0 million baryonic-|-DM par- 
ticles with mgas=m*= 2.5x10* and mDM=1.8xlO'' h~^MQ, 
egas=e*=2.8 and eDM=5.4 /i~^kpc. As a resolution test, the 
Virgo cluster was also simulated at this (eight times lower) 
mass and (two times lower) force resolution. This simulation 
will be denoted "VirgoLR". 

For all simulations, gravitational softening lengths were 
fixed in co-moving coordinates till z=6, subsequently in 
physical coordinates. 

For the simulation presented in this paper /wind=0.8, 
and an Arimoto-Yoshii IMF (of slope x—-l, shallower than 
the Salpeter slope x=-1.35) with mass limits [0.1-100] Mq 
was adopted. AGN driven winds were not invoked. More- 
over, the thermal conductivity was set to zero assuming that 
thermal conducti on in the ICM is highly suppressed by mag- 
netic fields (e.g.. lEttori fc Fabian! [2OO0I I . We note, that in 
relation to all relevant quantities presented in this paper, 
no significant difference is found between simulations with 
zero thermal conductivity and simulati ons invoking therm al 
conduction at 1/3 of the Spitzer level (| Romeo et ahllioosl ). 
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Table 1. Numerical and physical characteristics of the cluster simulations: mass of DM/gas/star particles and the respective gravitational 
(spline) softening lenghts; total number of particles and initial redshift of each run; virial mass, radius and temperature of simulated 
clusters (last three columns, referring to 2=0). 



run 










^gas 

[kpc/h] 




Ntot 




[1014 Mq] 


Rvir 

[Mpc] 


kT 
[keV] 


Coma 


180 


25 


25 


5.4 


2.8 


2.8 


950000 


19 


12.4 


2.9 


6 


Virgo 


23 


3.1 


3.1 


2.7 


1.4 


1.4 


2235000 


39 


2.8 


1.7 


3 


VirgoLR 


180 


25 


25 


5.4 


2.8 


2.8 


260000 


19 


2.8 


1.7 


3 


ComaDM 


26 






2.8 






4000000 


39 


12.4 


2.9 




VirgoDM 


26 






2.8 






1200000 


39 


2.8 


1.7 





2.1 Correcting for the effects of late 
quasi-stationary cooling flows 

After a period of major merging at z > 2 strong, quasi- 
stationary cooling flows develop at the centers of the clus- 
ters despite the strong, super-nova driven energy feedback 
to the IGM/ICM through galactic super- winds and the use 
of a fairly top-heavy IMF. As a consequence, stars continue 
to form steadily at the centers of the cDs (r <10 kpc) at 
rates of ~100 and ~500 Ma/ yr at ^=0 for Virgo and Coma , 
respectively. As discussed in ISommer- Larsen et al] (|2005h . 
this results in the cDs becoming too blue over time, with 
central B-R colours of ~1.0 at 2=0, rather than 1.4-1.5 ob- 
served for real cDs (note though that some star-formation is 
observed in many cDs at the base of the cooling flow — e.g., 
McNamara 2004). Given the significant gas cool-out rates, 
the increase in the stellar masses of the cDs since z ~2 is 
considerable: For the Virgo cD M«,i2=0.42, 1.0 and 1.4 in 
units of IO^^Mq, at z=2,l,0, respectively (inside of 20 kpc). 
For Coma, the corresponding masses are 0.95, 1.9 and 4.2. 
In real cDs, the late time cooling flows will be counteracted, 
most likely by the energy release related to the accretion of 
mass onto super-m assive black holes (e.g. jBower et al ]|2006l : 
ICroton et aLll2006l ). This will lead to a strong suppression 
of the late time star formation rate, such that the simulated 
cDs become too massive. As will be shown in the next sec- 
tion, the gravitational force of this excess mass influences 
the central dark matter profile, and must be corrected for. 

Assuming that the excess mass has been added adiabat- 
ically, i.e. on a timescale large compared to the stellar orbital 
timescale in the cD, the "correct" dark matter distribution 
can be recovered by slowly removing the excess stellar mass. 
This is accomplished as follows: 

In order to enable a direct comparison to the obser- 
vational results obtained for Abell 1703 we focus on dark 
matter profiles at 2=0.28, corresponding to t=10.2 Gyr. At 
t=8.2 Gyr, stars at r < 20 kpc from the center of the cD 
are identified. For each of the two simulated clusters two ad- 
ditional simulations are carried out — in one the masses of 
all identified stars formed since z=l are then gradually re- 
duced to zero over a period of 1 Gyr, in the other the masses 
of all identified stars formed since z=2. At the same time 
radiative cooling as well as star formation is switched off. 
Subsequently the additional simulations are then continued 
for 2 Gyrs more, still with radiative cooling and star forma- 
tion switched off. We refer to these additional simulations 
as Rzl and Rz2 respectively. 

The purpose of these additional simulations is to, with- 
out resorting to completely new simulations with AGN feed- 



back etc. included, determine the cluster dark matter profiles 
in the (more realistic) case where effects of late time cooling 
flows are counteracted since redshifts of either 1 or 2. 

2.2 Dark matter only simulations 

To further compliment the existing simulations and deter- 
mine the role of the inclusion of baryonic physics, we carried 
out dark matter only simulations of the two clusters using 
the same initial conditions as previously, but not splitting 
the particles into SPH and DM particles. 

The Virgo and Coma cluster were re-simulated using 1.2 
and 4.0 million DM particles, respectively. These simulations 
will be denoted "VirgoDM" and "ComaDM" in the follow- 
ing. The high resolution DM particles had mDM=2.6xlO* 
/i^iMq and eDM=2.8 h^^'kpc. 



Numerical and physical parameters for the simulations and 
the z=0 clusters are summarized in Table. [1] 



3 RESULTS 

3.1 Inner dark matter halo profiles 

To enable a precise determination of the galaxy cluster dark 
matter density profiles at 2=0.28, 21 time frames, with 
a time spacing of 0.1 Gyr (large enough to enable quasi- 
random sampling of the central halo profile) and centered 
on t=10.2 Gyr, were co-added. The DM profile power law 
index. 



where r is the cluster-centric radial distance, was then de- 
termined by averaging the co-added DM distribution over 
spherical shells. 

In Fig. [l] is shown by dotted lines the resulting a{r) 
(strictly speaking -a{r), but we shall neglect this distinc- 
tion in the following) for the dark matter only simulations. 
Moreover, a(r) is shown by thick solid curves for the orig- 
inal simulations including baryons, i.e. simulations includ- 
ing baryons, but without any corrections for late cooling- 
flows and related star formation. Finally, shown by thin 
solid lines, is a{r) for generalized NFW profiles (eq.[l]) with 
Qjv Fiy=a(0) = 1.09 and .92, as determined for A bell 1703 
bv iLimousin et all (j2008h and iRichard et all (j2009l l. respec- 
tively, and concentration parameter c=6. It is seen from the 



Pinching of dark matter haloes 5 




Figure 1. Logarithmic density slopes of the spherically averaged 
cluster dark matter distributions at 2:=0.28. Solid red curve shows 
result for the "Virgo" cluster including the cD uncorrected for ef- 
fects of late time cooling flows. Dotted red curve shows the result 
of the Virgo DM only simulation. The two red dashed lines show 
the results of the Rzl (lower) and Rz2 (upper) re-simulations, 
where effects of spurious late time cooling flows have been cor- 
rected for. The blue curves show the same for the "Coma" cluster. 
The two black solid curves show generalized NFW models (eq. [1]) 
with a]vFVl'=0-92 (upper) and 1.09 (lower), and concentration 
parameter c=r20o/rs=G. 

figure that a) all the simulations match the NFW profiles in 
the intermediate and outer parts of the haloes, r/r2oo ^ 0.07, 
fairly well, b) the dark matter only simulation of Virgo is 
matched quite well by an a(0) ~ 1.1 NFW profile and that 
of Coma by an q(0) ~ 0.9 profile, and c) the simulations 
with baryons result in too large values of a{r) (~ 1.5 — 1.6) 
in the inner parts, r/r2oo ^ 0.04 (the results are afi'ected by 
effects of gravitational softening at r/r2oo 0.007-0.008 — 
see section 3.2.1. for more detail). In order to assess whether 
this is due to baryonic "over" -pinching, i.e. due to the sim- 
ulated cDs being unrealistically massive, we compare the 
simulated and observed Virgo and Coma clusters: 

At z=0 the cDs in the simulated clusters have masses 
(inside of r=50 kpc — see below) of 7.4x10^^ and 
2.2xlO^^M0, for Coma and Virgo, respectively (at z=0.28 
the corresponding numbers are 5.1x10^^ and l.Sx 10^^ Af©). 
The real Coma and Virgo clusters both contain two domi- 
nant galaxies: NGC 4874 and NGC 4889 in Coma, and M86 
and M87 in Virgo. The two dominant galaxies in each clus- 
ter will eventually merge into one cD in each. The B-band 
luminosities of NGC 4874 ar id NGC 4889 are 1-5 xlO" and 
1.8xl0"Lo,s, respectively (jMichard fc Andreonll2008h . As- 
suming a B-band M/L of about eight (as found for the cDs in 
the present simulations based on the Arimoto-Yoshii IMF; 
the more standard Salpeter IMF results in similar M/L) 
the combined stellar mass of the two galaxies (at ^ ~ 0) 
is 2.6xlO^'^M0. Moreover, the simulated Coma cluster is 
a Tx ~6 keV cluster, whereas the real Coma cluster is a 



Figure 2. For re-simulations Rzl (solid curves) and Rz2 (dashed 
curves), at 2=0, are shown the cD cumulative stellar masses 
(Virgo: red curves; Coma: blue curves). 



Tx ~8 keV cluster. Assuming a standard scaling of mass 
T^''^ the cD in a 6 keV cluster should have a stellar mass 
of 1.7xlO^^M0. The B-band luminosities of M86 and M87 
are 4 .2x10^° and 6.6x 1O^°L0,b, respectively (|Bingelli et all 
Il985h . The stellar masses of the two galaxies thus add up to 
about 8.6x10"Mq. 

The simulated cDs are hence about 3-4 times more mas- 
sive than observed, indicating that the simulated DM haloes 
are too contracted compared to real ones, which at least 
qualitatively could explain some of the differences seen in 

Fig.m 

To understand the differences seen in Fig.[T]more quan- 
titatively we now turn to the re-simulations aimed at cor- 
recting for the effects of the late cooling fiows. In Fig. [T] 
we show by thick dashed lines a{r) for the simulations with 
baryons, where the mass of all stars in the central cD, formed 
since either Zf=l or 2, has been adiabatically reduced to zero 
— we shall denote these two sets of re-simulations by Rzl 
and Rz2 in the following. For the inner parts of Coma, a 
large effect on a(r) results already from removing central 
stars formed since zf = l. This is related to the fact that the 
mass of the cD more than doubles since z=l in the original 
simulation including baryons. For Virgo a strong effect is 
seen only between the original simulation and Zf=2 correc- 
tion. This is also reasonable, since the mass of the cD more 
than doubles going from z—2 to 1, but then only increases 
by additional 40% going to redshift 0. In order to assess 
which of the Rzl and Rz2 re-simulations, if any, result in 
cDs which are realistic, we now derive the properties of the 
resulting cDs. 
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3.2 Properties of the simulated cD galaxies 

3.2.1 cD stellar masses 

In Fig. |2]is shown the cumulative stellar mass for the Coma 
and Virgo cDs from the z/ = l and 2 adiabatically corrected 
simulations, respectively. It is seen, that at r—50 kpc, cor- 
responding to approximately the 25 B-mag/arcsec'^ (opti- 
cal) isophote (see below), for Zf = l the cumulative stellar 
masses are 3.1 and I.SxIO^^Mq for Coma and Virgo, re- 
spectively. For Zf=2 the corresponding numbers are 2.1 and 
0.74x lO^'^M©. Comparison to the observational estimates 
above, indicates that the z/=2 cDs provide a better match 
to reality than the z/ = l ones. 

An alternative c omparison can be ma de to the cD of 
Abell 1703, for which [Limousin et all (|2008l ) estimate a stel- 
lar mass of 1.3±0.3 Mq inside of r=30 kpc. The tempera- 
ture of the Abell 1703 ICM is uncertain, but (also uncer- 
tain) estimates of the virial mass of Abell 1703 indicate that 
the temperature is likely at least as large as for the simu- 
lated Coma cluster (Riemer-S0rensen, private communica- 
tion, Oguri et al 2009). For the Coma Rzl and Rz2 cDs, at 
2=0.28, we find stellar masses inside of r=30 kpc of 2.2 and 
1.2xlO^^M0. A gain, the stellar mass appears realistic for 
the Rz2 Coma cD, but too large for the Rzl cD. 




Figure 3. Stellar oxygen abundanee as a funetion of projected 
radial distance from the center of the cD for re-simulations Rzl 
and Rz2, at z=0. The curves are labeled as in Fig. |2] 



3.2.2 cD effective radii 

Next, we determine the effective radii for the cDs at z—Q. To 
this end, the projected B-band surface brightness profiles are 
determined for the cDs (see Romeo et al. 2005 for details on 
calculating broad band photometric properties of the star 
particles): First, a B-band surface brightness map is con- 
structed by projecting the star-particle distribution within 
a cube of side-length L=200 kpc, centered on the cD, along 
the three cardinal directions and then superposing the three 
maps (larger values of L result in similar results). Second, a 
radial surface profile is constructed by azimuthally averaging 
over the combined image. Third, the radial profile obtained 
is divided by a factor of three. At R ~50 kpc the B-band 
surface brightness for both cDs drops to 25 mag/arcsec'^ 
ISommer-Larsen et al.l [20051 ). which is taken to be the ra- 
dial extent of the cD itself. -Rofi is then determined as the 
projected radius enclosing half of the cD B-band luminos- 
ity. For the Coma and Virgo cDs in the Rzl simulations, 
Ref[=8 and 7 kpc is found. The corresponding values for the 
Rz2 simulations are _Roff = 16 and 11 kpc. For a large (ob- 
servat ional) sample of SDSS galaxies (at z ^0). IShen et al.l 
l|2003h find for early type galaxies of stellar masses 3.1 and 
1.3xlO^^M0 median effective radii of ~25 and 16 kpc, and 
for stellar masses of 2.1 and 0.74x 10^^ Af©, Rc« ~20 and 
13 kpc. It follows that the Rzl cDs have too small effective 
radii for their mass, whereas the Rz2 cDs have reasonable 
effective radii compared to observations. 



3.2.3 cD metal abundances and alpha/Fe ratios 

In Fig.|3]is shown, for cDs in the Rzl and Rz2 re-simulations, 
the projected average stellar oxygen abundance as a func- 
tion of projected radius R. In these alpha-element enhanced 
systems (see below), [0/H] can be taken as a measure of 



the total abundance [Z/H]. It is seen that the oxygen abun- 
dance is super-solar all the way to 50 kpc projected radius. 
For the Rzl simulations, the integrated (global) cD abun- 
dances are [O/H]=0.59 for both clusters. For the Rz2 simu- 
lations, the integrated abundances are [O/H] =0.40 and 0.26 
for Co ma and Virgo, respectively. Following [Thomas et al.l 
(|2005h . observed early type galaxies of stellar masses 3.1 
and 1.3xlO^^M0 have median <[Z/H]>=0.38 and 0.33, re- 
spectively, and for stellar masses of 2.1 and 0.74x10"Mq, 
<[Z/H]>=0.36 and 0.30. So the Rzl cDs have too large 
abundances compared to what is observed, whereas the Rz2 
cDs provide a good match to observations. Taking [O/Fe] 
as a proxy for alpha/Fe, the results for the Rzl and Rz2 
simulations are similar, [O/Fe] ~0. 3-0. 35. These values are 
within the range o bserved for large early type galaxies 
([Thomas et al.ll2005l ). 

3. 2.4 cD U-V colours 

For 2 ~ early type galaxies, the rest-frame U-V broadband 
colour allows sampling of the 4000 A break, and therefore 
is particularly sensitive to age and metallicity variations of 
the stellar populations. In order to check whether the simu- 
lated cDs fall on the observational "red sequence" for early 
type galaxies, we hence calculate the U-V colours for the 
cDs. In Fig. |4] is shown, for cDs in the Rzl and Rz2 re- 
simulations, the U-V colour as a function of projected ra- 
dius R. In the inner parts of the galaxies, all the cDs dis- 
play a negative colour gradient, qualitatively consistent with 
what is observed for large early type galaxies. As the mean 
stellar age is approximately constant with projected radius 
R, the reason for the gradient in colour is the metallicity 
gradient seen in Fig. |31 This is also the generally accepted 
explanation f or the observed colo ur gradients in early type 
systems {e.g. [Tamura et al.|[2000l l. For the Rzl cDs, inte- 




Figure 4. U-V colour as a function of projected radial distance 
from the center of the cD for re-simulations Rzl and Rz2, at z=0. 
The curves are labeled as in Fig. |2] 



grated (global) U-V colours of 1.31 and 1.27 are found for 
Coma and Virgo, respectively. For the R2z cDs, the corre- 
sponding numbers are U-V=1.29 and 1.28. The two Rzl cDs 
have My=-24.6 and -23.7 for Coma and Virgo, respectively. 
For the Rz2 cDs, the corresponding numbers are -23.8 and 
-23.0 . Comparing to the observational U-V vs. Mv sequence 
{e.g. iRomeo et al.|[200^ ). early type galaxies of such abso- 
lute V magnitudes have a median U-V~1.4. We hence find 
that the Rzl and Rz2 cDs have very similar U-V colours, 
and that these are slightly bluer than the observed average. 
This is, at least partly, due to the large, somewhat bluer 
stellar envelope surrounding the cDs, cf. Fig. |4l 



3.2.5 cD central physical stellar densities 

From the results presented above it is possible to estimate 
the stellar core density of the cDs, defined as 



Pc,* 



4^/3 



(3) 



For the Rzl cDs we find pc,*=6 and 4x10* Mo/kpc^ 
for Coma and Virgo, respectively. For the Rz2 cDs, the 
corresponding numbers are 5 and 4x10^ Mo/kpc^. Ex- 
trapolating the ob s ervat ional estimates of Pc,*, given in 
Ivan Dokkum et "aP l)2008h on the basis of SDSS data, to 
galaxy stellar masses of ~ 10^^ Mq, the Rz2 cD core den- 
sities appear the most reasonable, although no strong con- 
clusions can be made on the basis of the available data. 



3.3 How much are galaxy cluster dark matter 
haloes pinched by the central cDs? 

Based on the results given in the previous subsection, it is 
clear that the cDs formed in the Rz2 re-simulations give a 
much better match to reality than the ones formed in the 



Figure 5. Logarithmic density slopes for the Rz2 re-simulations. 
Results for total matter densities are shown by thick, short-dashed 
linos (red: Virgo; blue: Coma). The thin long-dashed lines shows 
the result for the dark matter (as in Fig.[T]l. The two black solid 
curves show generalized NFW models (cq.[l]) with qjvfw=1-0 
and 1.3, and concentration parameters c=r20o/''s=4 and 6, re- 
spectively. 



Rzl simulations, and in fact, in general, meet all observa- 
tional constraints. We shall hence in the following focus on 
the Rz2 re-simulations, and the comparison of these to the 
dark matter only simulations. 



3.3.1 Resolution limitations due to gravity softening 

The Coma dark matter only simulation has eDM=2.8 
/i~^kpc, and the gravity force of the central dark matter par- 
ticles is hence purely Newtonian (i.e. un-softened) at r ~8 
kpc. At 2=0.28, r2oo for the Coma cluster is 1820 kpc, which 
means that the simulation result is affected by resolution ef- 
fects inside of r/r2oo —0.004. The Coma simulation with 
baryons has eDM=5.4 /i~^kpc, so although the star and gas 
particle softening lengths are considerably smaller, the re- 
sulting dark matter profile may be affected by gravity soft- 
ening inside of r/r2oo —0.008. The Virgo simulations with 
and without baryons have eDM=2.7 and 2.8 /i~^kpc, respec- 
tively. As r2oo for this cluster at z=0.28 is 1105 kpc, the 
simulation results may be affected by gravity softening in- 
side of r/r2oo —0.007. Finally, the VirgoLR simulations have 
eDM=5.4 and 5.6 /i~^kpc, respectively, resulting in an inner 
radial resolution limit of r/r2oo —0.014. 

In the following we shall only discuss simulation results 
outside of the above gravity softening resolution limits. 



3.3.2 Steepening of the inner dark matter profiles 

As can be seen from Fig. [T] the dark matter density pro- 
file obtained in the Coma dark matter only simulation, can 
be fairly well described by a generalized NFW profile of 




Figure 6. For the Rz2 rc-simulations is shown a(r) for the da,rh 
matter alone (thick short-dashed lines), and for dark matter and 
ICM gas combined (thin long-dashed lines). Results for Coma 
are shown by blue curves, for Virgo by red curves. The two black 
curves are the NFW models shown in Fig. [T] 



Figure 7. Comparison of the Virgo and VirgoLR simulations, at 
2=0. The red curves display the results for Virgo, labeled in the 
same way as in Fig. [T] The black curves show the correspond- 
ing results for VirgoLR. The statistical uncertainty for the Vir- 
goLR data is Ao ~0.05; for Virgo about a factor of \/8 less. The 
light blue curve shows a generalized NFW models (eq.[l]) with 
ojVFW=l-09 and concentration parameter c=r2oo/^s=S. 



aNFW ~ a(0)~0.9. Comparing it to what is obtained for 
the Rz2 re-simulation of Coma the steepening of the in- 
ner dark matter profile, caused by the cD, is quite mod- 
erate, Aa ~0.1. For the Virgo cluster, the dark matter 
only profile is also fairly well described by a generalized 
NFW profile, this time of aNFW = a(0)~l.l. Again, the 
Rz2 re-simulation only leads to a very moderate steepening, 
Aq ~0.1, relative to the dark matter only simulation. So 
for reahstic cDs, the steepening of the dark matter profile in 
galaxy clusters is very moderate, Aq ~0.1. This is the most 
important result presented in this work. 



3.3.4 Dark matter + gas density profiles 

In performing the dark matter vs. baryonic matter decompo- 
sition mentioned above, one typically assumes that the bary- 
onic component is dominated by the stellar mass of the cD, 
for which the mass distribution can be fairly well modeled. 
It is hence important to show, that neglecting the hot gas 
(ICM) component, for which the mass distribution is much 
more uncertain, does not significantly bias the estimate of 
the DM distribution. In Fig. [S] is shown, for the Rz2 simu- 
lations, a(r) for the dark matter alone, as well as for dark 
matter and hot ICM gas combined. As can be seen, the effect 
on a{r) of including the ICM is very modest, Aq <0.05. 



3.3.3 Total matter density profiles 

As gravitational lensing probes the total matter (surface) 
density, in order to deduce dark matter profiles observation- 
ally, one has to do a dark matter vs. baryonic matter decom- 
position. This obviously adds some uncertainty (statistical 
as well as possibly systematic) to the estimate of the dark 
matter density profile. To bypass this, it is clearly also of 
interest to calculate total matter density profiles from the 
simulations. These can then more directly be compared to 
observational finding, based on lensing. In Fig. |5] is shown 
the result for the Rz2 simulations. Both the Virgo and the 
Coma profile are characterized by a rather sharp break at 
r/r2oo —0.05-0.06. Inside of this the profiles are rather fiat, 
with a{r) ~ Q-NFW;DMoniy + 0.5. The profiles are not well 
described by generalized NFW models, as exemplified in the 
figure. 



3.4 A numerical resolution test 

In Fig. [7] is shown a{r) for the sets of Virgo simulations 
as well as VirgoLR simulations. Statistical uncertainties on 
the VirgoLR results are Aq ~0.05, and on the Virgo results 
about a factor \/8 less. Note that the VirgoLR results are 
affected by gravity softening at about r/r2oo —0.014. Out- 
side of this region, there is reasonably good agreement be- 
tween the two sets of simulations, indicating that the results 
presented in Fig. [1] are not seriously affected by numerical 
resolution limitations outside of r-/r2oo —0.007-0.008. 



4 DISCUSSION 

Fitting the dark matter density profile of Abell 1703 by 
a generalized NFW profile, different authors find q(0) be- 
tween 0.8 and 1.2 (|Limousin et all I2OO8I : [Saha et all l2009l : 
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iRichard et "aill2009l : lOguri et all 120091 '). Although we in this 
paper are analyzing simulations of just two clusters, and al- 
though other observational values of Q:(0) have previously 
been reported in the literature, it is obviously interesting 
that we for the two clusters simulated, including the effects 
of the central cD on the dark matter halo, find a{0) ~1.0- 
1.2. 

Given that LCDM haloes undoubtedly have a range 
of inner slopes, and given that we are probing only two 
haloes, an even more interesting result is that for both clus- 
ters, the steepening of the inner slopes is both moderate 
and consistent, Aa ^0 .1. For Milky Way sized galaxies, as 
well as smaller galaxies, iGustafsson et aL (|2006l l found much 
larger effects of baryonic pinching of the dark matter haloes, 
Aq ~0.6. To understand the reason for this difference, we 
show in Fig. |8] for the Rz2 re-simulations at 2=0.28, the 
cumulative masses of dark matter, stars and stars-|-gas. As 
can be seen, the cumulative masses of baryons and dark 
matter are equal at req —7 and 9 kpc, for Virgo and Coma 
respectively. In both cases this is well inside of the effective 
radius , i.e. the charact e ristic scale of the galaxy. In con- 
trast, IGustafsson et af] (|2006l ) found req ~10 kpc for the 
individual (field) galaxies considered. These galaxies have 
characteristic scales of ~2-4 kpc, so the chief difference is 
that the cD galaxies considered in this paper are much less 
baryon dominated than normal (substantially smaller) field 
galaxies. In addition, one should note that the present sim- 
ulations can only probe the increase in dark matter inner 
profile steepening outside of ~8 and 16 kpc, for Virgo and 
Coma respectively (whereas the innermost, baryon domi- 
nated regions are about twice as well resolved due to the 
shorter gravitational softening lengths of the star and gas 
particles) . 

It is also worth noting that lCnedin et all (|2004l l found 
a consideably larger steepning of galaxy cluster inner DM 
profiles due to baryonic pinching, on the basis of their sim- 
ulations invoking baryonic physics. This, however (as ac- 
knowledged by these authors), is due to the spurious effects 
of "over-cooling" , discussed and corrected for in the present 
work. 

An underlying assumption in the Rzl and Rz2 re- 
simulations is that the rate of change of cD gravitational po- 
tential, caused by stars either being formed m situ or being 
added through accretion since Zf — 1 or 2, is slow (adiabatic), 
compared to the dynamical rate. Although this is true for 
the stars forming m situ at the base of the cooling flows, 
in the Coma simulation ~1:3 and 1:4 mergers take place at 
z ~1.4 and z ~0.7, respectively, and in the Virgo simulation 
~1:4 mergers take place, at z ~1.5 and 1.0, respectively. In 
such mergers, the adiabatic approximation is not fulfilled, 
and more realistic simulations, including effects of super- 
massive black hole growth and AGN feedback, would have 
to be undertaken to check whether the results we find from 
the adiabatic stellar mass removal re-simulations are in fact 
fully correct. However, the fact that the Rz2 cDs appear to 
be realistic in terms of stellar mass, linear extent, metal- 
licity, etc. can be considered very strong evidence that the 
inner dark matter profile steepening found for the Rz2 simu- 
lations are entirely reasonable. It is obvious that the primary 
drivers of the baryonic pinching are a) cD stellar mass, and, 
to a lesser extent, b) the cD linear extent. 
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Figure 8. For the Rz2 runs, at 2=0.28, arc shown cumulative 
mass of dark matter (solid curves; Coma: thick, Virgo: thin), stars 
(dotted curves; Coma: blue, Virgo: red) and stars-|-gas (dashed 
curves; Coma: blue, Virgo: red). 



5 CONCLUSION AND OUTLOOK 

We have in this paper presented results of cosmological 
(ACDM) TreeSPH simulations of the formation and evolu- 
tion of two galaxy clusters. In order to quantify the effect of 
baryonic "pinching" of galaxy cluster dark matter haloes, we 
have simulated the same clusters with and without baryons 
included. 

The hydrodynamical simulations include a number of 
important physical processes, such as star formation, chem- 
ical evolution with non-instantaneous recycling, metallicity 
dependent radiative cooling, strong star-burst, driven galac- 
tic super-winds and the effects of a meta-galactic UV field, 
including simplified radiative transfer. The two clusters have 
T~3 and 6 keV, respectively, and both host at 2: ~ a 
prominent, central cD galaxy. 

Comparing the simulations without and with baryons, 
it is found that for the latter, the inner DM density profiles, 
r < 50-100 kpc, steepen considerably: Aa ~0.5-0.6, where 
-Q is the logarithmic DM density gradient. This is primar- 
ily due to the central stellar cDs becoming very massive, as 
a consequence of the onset of late time cooling fiows and 
related star formation. Once these spurious cooling fiows 
have been corrected for, and the cluster gravitational po- 
tentials dynamically adjusted, much smaller pinching effects 
are found: Aa ~0.1. 

The main result of the paper is hence that our work 
strongly supports the notion that the steepening of the inner 
dark matter density profile, due to baryonic pinching, is very 
moderate. This finding is at stark variance with findings for 
individual field galaxies. The main reason for this is that 
galaxy clusters, even in the presence of central dominant cD 
galaxies, are dark matter dominated outside of ~7-8 kpc. As 
the effective radii of the two cDs considered are 11 and 16 
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kpc, this result indicates that cD galaxies are dark matter 
dominated at the effective radius. 

Including the effects of baryonic pinching, central slopes 
of a ~1.0 and 1.2 are found for the DM in the two clus- 
ters, interestingly close to recent observational findings for 
Abell 1703 based on gravitational lensing studies. 

For the simulations with baryons, the inner density 
profile of DM+ICM combined is found to be only very 
marginally steeper than that of the DM. The total matter 
inner density profiles, however, are found to be Aa ~ 0.5 
steeper than the inner profiles in the dark matter only sim- 
ulations. The total matter density profiles are not well de- 
scribed by generalized NFW models. 

Given that we have presented results for just two sim- 
ulated clusters in this paper, it is clearly of importance to 
increase the sample size. On the observational side, strong 
lensing constraints for more uni-modal clusters should be 
obtained. We plan to present such results in forthcoming 
papers. 
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